body .main-container {
max-width: 100% !important;
width: 100% !important;
}
body {
max-width: 100% !important;
}
body, td {
font-size: 16px;
}
code.r{
font-size: 14px;
}
pre {
font-size: 14px
}Let us set some global options for all code chunks in this document.
knitr::opts_chunk$set(
message = FALSE, # Disable messages printed by R code chunks
warning = FALSE, # Disable warnings printed by R code chunks
echo = TRUE, # Show R code within code chunks in output
include = TRUE, # Include both R code and its results in output
eval = TRUE, # Evaluate R code chunks
cache = FALSE, # Enable caching of R code chunks for faster rendering
fig.align = "center",
out.width = "100%",
retina = 2,
error = TRUE
)
rm(list = ls())
set.seed(1982)Let us now load some required libraries.
# Load required libraries
# inla.upgrade(testing = TRUE)
# remotes::install_github("inlabru-org/inlabru", ref = "devel")
# remotes::install_github("davidbolin/rspde", ref = "devel")
# remotes::install_github("davidbolin/metricgraph", ref = "devel")
library(INLA)
library(inlabru)
library(rSPDE)
library(MetricGraph)
library(plotly)
library(dplyr)
library(tidyr)
library(MASS)
library(glmnet)
library(car)
library(ggplot2)
library(plotly)
library(sf)
library(here)We now define some R functions with self-explanatory
names.
# Functions to manipulate the data
standardize <- function(vector) {return((vector - mean(vector)) / sd(vector))}
minmax_scaling <- function(x) {return((x - min(x)) / (max(x) - min(x)))}
convert_to_binary <- function(input_vector) {return(ifelse(input_vector != 0, 1, 0))}The following chunk uses the function
creates_covariates_on_mesh() defined on
Preprocessing/28Window.case.file4.addcovariatesonmesh.R to
build a mesh (on the same graph as the one used in this file) and
compute the value of all available covariates on the mesh nodes.
Go to Preprocessing 8.0.1 Raw
code to see the body of the
creates_covariates_on_mesh() function.
h = 0.05
# source(here("Preprocessing/28Window.case.file4.addcovariatesonmesh.R"))
# creates_covariates_on_mesh(h)
# Sys.sleep(5)
################################################################################
################################# PREPARE THE DATA #############################
################################################################################
# loading the data
load(here("Data_files/data_on_graph_with_covariates_no_consecutive_zeros_partialtomtom.RData"))
load(here("Graph_objects/graph_construction_28_03_2024partialtomtomwhichlonglatsf.RData"))
load(here("Data_files/data_on_mesh_with_covariates_partialtomtom.RData"))
data = data_on_graph_with_covariates %>%
mutate(across(starts_with(c("class_", "upto")), list(ind = convert_to_binary))) %>%
mutate(across(c("bus", "signal", "stop", "crossing"), ~round(., 5))) %>%
mutate(across(c("SpeedLimit", "density_per_hour"), standardize)) %>% dplyr::select(-datetime)
mesh = data_on_mesh_with_covariates %>%
mutate(across(starts_with(c("class_", "upto")), list(ind = convert_to_binary))) %>% # this creates new columns
mutate(across(c("bus", "signal", "stop", "crossing"), ~round(., 5))) %>%
mutate(across(c("SpeedLimit", "density_per_hour"), standardize))Let us add the speed observations to the graph and build the mesh.
Observe that mesh already has the value of the available
covariates on the mesh nodes.
sf_graph$add_observations(data = data, group = "day", normalized = TRUE, clear_obs = TRUE)
sf_graph$build_mesh(h = h)
Sys.sleep(5)stat.time.ini <- Sys.time()
################################################################################
################################# STATIONARY MODEL #############################
################################################################################
rspde_model_stat <- rspde.metric_graph(sf_graph,
parameterization = "matern",
nu = 0.5)
data_rspde_bru_stat <- graph_data_rspde(rspde_model_stat,
repl = ".all",
name = "field")
st.dat.rep.stat <- inla.stack(
data = data_rspde_bru_stat[["data"]]["speed"],
A = list(data_rspde_bru_stat[["basis"]],1),
tag = "est",
effects =
list(c(data_rspde_bru_stat[["index"]], list(Intercept = 1)),
list(SpeedLimit = data_rspde_bru_stat[["data"]]["SpeedLimit"])
)
)
f_stat = speed ~ -1 + Intercept + SpeedLimit + f(field, model = rspde_model_stat, replicate = field.repl)
rspde_fit_stat <-
inla(f_stat,
data = inla.stack.data(st.dat.rep.stat),
family = "gaussian",
control.predictor = list(A = inla.stack.A(st.dat.rep.stat))
)
stat.time.fin <- Sys.time()
print(stat.time.fin - stat.time.ini)## Time difference of 15.23668 secssummary(rspde_fit_stat)##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " keep = keep, working.directory = working.directory,
## silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
## debug, .parent.frame = .parent.frame)" )
## Time used:
## Pre = 0.446, Running = 5.8, Post = 0.98, Total = 7.22
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## Intercept 20.883 0.374 20.163 20.877 21.638 20.877 0
## SpeedLimit 3.060 0.280 2.517 3.061 3.597 3.062 0
##
## Random effects:
## Name Model
## field CGeneric
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.012 0.000 0.012 0.012
## Theta1 for field 2.853 0.070 2.717 2.852
## Theta2 for field -1.089 0.217 -1.508 -1.091
## 0.975quant mode
## Precision for the Gaussian observations 0.012 0.012
## Theta1 for field 2.993 2.849
## Theta2 for field -0.656 -1.102
##
## Marginal log-Likelihood: -55592.19
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')fit.rspde = rspde.result(rspde_fit_stat, "field", rspde_model_stat)
summary(fit.rspde)nonstat.time.ini <- Sys.time()
################################################################################
############################# NON STATIONARY MODEL #############################
################################################################################
B.sigma = cbind(0, 1, 0, mesh$SpeedLimit, 0)
B.range = cbind(0, 0, 1, 0, mesh$SpeedLimit)
rspde_model_nonstat <- rspde.metric_graph(sf_graph,
start.theta = c(fit.rspde$summary.log.std.dev$mode,
fit.rspde$summary.log.range$mode,
rep(0, 2)),
B.sigma = B.sigma,
B.range = B.range,
parameterization = "matern",
nu = 0.5)
data_rspde_bru_nonstat <- graph_data_rspde(rspde_model_nonstat,
repl = ".all",
name = "field")
st.dat.rep.nonstat <- inla.stack(
data = data_rspde_bru_nonstat[["data"]]["speed"],
A = list(data_rspde_bru_nonstat[["basis"]],1),
tag = "est",
effects =
list(c(data_rspde_bru_nonstat[["index"]], list(Intercept = 1)),
list(SpeedLimit = data_rspde_bru_nonstat[["data"]]["SpeedLimit"])
)
)
f_nonstat = speed ~ -1 + Intercept + SpeedLimit + f(field, model = rspde_model_nonstat, replicate = field.repl)
rspde_fit_nonstat <-
inla(f_nonstat,
data = inla.stack.data(st.dat.rep.nonstat),
family = "gaussian",
control.predictor = list(A = inla.stack.A(st.dat.rep.nonstat))
)
nonstat.time.fin <- Sys.time()
print(nonstat.time.fin - nonstat.time.ini)## Time difference of 19.44462 secssummary(rspde_fit_nonstat)##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " keep = keep, working.directory = working.directory,
## silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
## debug, .parent.frame = .parent.frame)" )
## Time used:
## Pre = 0.2, Running = 12.4, Post = 1.02, Total = 13.6
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## Intercept 20.799 0.405 20.005 20.796 21.604 20.796 0
## SpeedLimit 2.863 0.273 2.332 2.862 3.404 2.863 0
##
## Random effects:
## Name Model
## field CGeneric
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.012 0.000 0.012 0.012
## Theta1 for field 3.011 0.114 2.784 3.012
## Theta2 for field -0.722 0.292 -1.302 -0.720
## Theta3 for field 0.371 0.119 0.133 0.372
## Theta4 for field 0.741 0.232 0.279 0.743
## 0.975quant mode
## Precision for the Gaussian observations 0.012 0.012
## Theta1 for field 3.234 3.015
## Theta2 for field -0.152 -0.713
## Theta3 for field 0.602 0.375
## Theta4 for field 1.193 0.751
##
## Marginal log-Likelihood: -55598.78
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')summary(rspde.result(rspde_fit_nonstat, "field", rspde_model_nonstat))save.image(here(paste0("Models_output/", rmarkdown::metadata$title, ".RData")))